rm(list = ls())

setwd('/path/to/replication/')

library(data.table)
library(ggplot2)
#library(devtools)
#devtools::install_github("naoki-egami/DIDdesign", dependencies = TRUE)
library(DIDdesign)

# Figure B9: Covid-19 death and unemployment effects on Asian hate crimes: Accounting for possible violations to the common trends assumption
load('./data/panel_month_dummies.RData')

# Figure B9 (a)
if (!file.exists('./output/doubledid_death_unbalanced_panel.RData')) {
  
  set.seed(848)
  fit_unbalanced_panel <- did(
    formula = hc_pc_asians ~ excess_deaths_jan_d + covid2,
    data    = panel,
    id_time = 'month',
    is_panel= FALSE,
    option  = list(n_boot = 200, parallel = TRUE, id_cluster = "panelvar", lead = c(0,1))
  )
  save(fit_unbalanced_panel, file = './output/doubledid_death_unbalanced_panel.RData')
}  

fit_unbalanced_panel <- get(load('./output/doubledid_death_unbalanced_panel.RData'))
results <- data.table(fit_unbalanced_panel$estimates)
results <- merge(results, data.table(lead = c(0,1), month = panel[covid2==1, unique(month)]), by = 'lead')
results[,c('conf.low', 'conf.high'):=.(estimate-1.96*std.error, estimate+1.96*std.error)]
results <- results[estimator!='sDID']

ggplot(results,aes(x = month, y = estimate, color = estimator, group = estimator)) +
  geom_point(position=position_jitterdodge(seed=196)) +
  geom_hline(yintercept = 0, colour='#999999', linetype='longdash') +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), position=position_jitterdodge(seed=196)) +
  scale_color_grey(end = 0.6) +
  scale_x_date(breaks = "1 month", minor_breaks = "1 month", date_labels = '%Y-%m',
               limits = c(as.Date('01-01-2020', format='%m-%d-%Y'), max = as.Date('04-01-2020', format='%m-%d-%Y')),
               expand=c(0,0)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position="bottom", legend.title = element_blank(), legend.text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"
        ),axis.title.y = element_text(size = rel(1.2)),axis.text.y=element_text(size=12)) +
  xlab(NULL) +
  ylab('Effect of health threat on Asian hate crimes')
ggsave('./robustness_double_did_death.pdf', width=6, height=4.85)


# Figure B9 (b)
if (!file.exists('./output/doubledid_unbalanced_panel.RData')) {
  
  set.seed(848)
  fit_unbalanced_panel <- did(
    formula = hc_pc_asians ~ losing_income_d + covid2,
    data    = panel,
    id_time = 'month',
    is_panel= FALSE,
    option  = list(n_boot = 200, parallel = TRUE, id_cluster = "panelvar", lead = c(0,1))
  )
  save(fit_unbalanced_panel, file = './output/doubledid_unbalanced_panel.RData')
}

fit_unbalanced_panel <- get(load('./output/doubledid_unbalanced_panel.RData'))

# plot results
results <- data.table(fit_unbalanced_panel$estimates)
results <- merge(results, data.table(lead = c(0,1), month = panel[covid2==1, unique(month)]), by = 'lead')
results[,c('conf.low', 'conf.high'):=.(estimate-1.96*std.error, estimate+1.96*std.error)]
results <- results[estimator!='sDID']

ggplot(results,aes(x = month, y = estimate, color = estimator, group = estimator)) +
  geom_point(position=position_jitterdodge(seed=196)) +
  geom_hline(yintercept = 0, colour='#999999', linetype='longdash') +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), position=position_jitterdodge(seed=196)) +
  scale_color_grey(end = 0.6) +
  scale_x_date(breaks = "1 month", minor_breaks = "1 month", date_labels = '%Y-%m',
               limits = c(as.Date('01-01-2020', format='%m-%d-%Y'), max = as.Date('04-01-2020', format='%m-%d-%Y')),
               expand=c(0,0)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position="bottom", legend.title = element_blank(), legend.text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"
        ),axis.title.y = element_text(size = rel(1.2)),axis.text.y=element_text(size=12)) +
  xlab(NULL) +
  ylab('Effect of unemployment threat on Asian hate crimes')
ggsave('./robustness_double_did.pdf', width=6, height=4.85)

